Matrix algorithm for solving Schrodinger equations with position-dependent mass or 

complex optical potentials 
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We represent low dimensional quantum mechanical Hamiltonians by moderately sized finite ma- 
trices that reproduce the lowest 0(10) boundstate energies and wave functions to machine precision. 
The method extends also to Hamiltonians that are neither Hermitian nor PT symmetric and thus 
allows to investigate whether or not the spectra in such cases are still real. Furthermore, the ap- 
proach is especially useful for problems in which a position-dependent mass is adopted, for example 
in effective-mass models in solid-state physics or in the approximate treatment of coupled nuclear 
motion in molecular physics or quantum chemistry. The performance of the algorithm is demon- 
strated by considering the inversion motion of different isotopes of ammonia molecules within a 
position-dependent mass model and some other examples of one- and two-dimensional Hamiltoni- 
ans that allow for the comparison to analytical or numerical results in the literature. 



I. INTRODUCTION 

In a number of cases complex systems can be treated 
in a simplified manner, if a position-dependent mass is 
introduced. One prominent example is the concept of 
an effective electron mass in solid-state physics where 
position-dependent masses provide a way to obtain cor- 
rections to the simplest approach in which a constant ef- 
fective mass is used (see, e. g., [l| and references therein). 
As a consequence, Schrodinger equations with a position- 
dependent mass have been considered in various contexts, 
for example to study electronic properties of semiconduc- 
tors 0, He clusters Q, or super-lattice band structures 
1- 

Also in molecular physics or theoretical chemistry 
position-dependent masses can occur, if high-dimensional 
nuclear motion is described using a lower dimensional ef- 
fective Hamiltonian. A rather well known example is the 
theoretical description of the inversion motion (umbrella 
mode) in ammonia molecules (NH3) in which the nitro- 
gen atom is moving from one side of the plane formed by 
the three hydrogen atoms to the other. In fact, this in- 
volves a collective motion, since the hydrogen-nitrogen 
bonds change their lengths while the nitrogen atom 
moves. This is thus an example of strongly coupled vibra- 
tional modes, in this case of the symmetric bending and 
stretching modes. While the main physics of this motion 
can be captured by an effective one-dimensional double- 
well potential, an improved approximation is achieved by 
introducing a position-dependent mass while still main- 
taining a one-dimensional treatment. The inversion mo- 
tion of ammonia and its one-dimensional model has been 
studied extensively in several papers, e.g. .in [B-Q, and 
has been exploited for the ammonia maser [9| . Especially 
Aquino et al. obtained very accurate results in com- 
parison to experiment by reducing the two-dimensional 
problem to a one-dimensional problem with a position- 
dependent reduced mass. Note that this ansatz is de- 
duced from classical arguments and thus leads to quan- 
tum mechanical operator ordering ambiguities, since it is 
not clear anymore how to order the operators in the ki- 



netic term p 2 /2m. The accurate and efficient solution of 
a Schrodinger equation with a position-dependent mass is 
non-trivial and various efforts were made even recently to 
find analytical or numerical solutions, e.g., in (ll. ITol. fill] . 

In this work a matrix method is introduced that de- 
termines efficiently and accurately the bound states of 
Schrodinger equations with a position-dependent mass. 
The algorithm is very flexible and allows to con- 
sider different symmetrized or non-symmetrized forms 
of the kinetic-energy operator. Furthermore, also non- 
Hermitian Hamiltonians can be treated, including com- 
plex ones. Non-Hcrmitian PT-symmetric Hamiltonians 
have recently stirred some interest [12], EH in connec- 
tion with the question under which circumstances the 
spectrum may still remain real. In addition to that, non- 
Hermitian Hamiltonians (also without PT-symmetry) oc- 
cur for example in the context of complex optical poten- 
tials (see, e.g., (ME!)- 

After a brief discussion of Hamiltonians with a 
position-dependent mass or PT symmetry in Sec. |H] 
the method is presented in Sec. Mil In Sec. IIVI the 
performance of the method is discussed for a number 
of examples. This includes a study of the inversion 
mode of NH 3 (Sec. IP7A")) and ND 3 (Sec. ITVBl) . the 
Morse potential (Sec. IIVD[) . a harmonic oscillator with 
position-dependent mass (Sec. IIVE[) . and examples of 
PT-symmetric and non-symmetric non- Hermitian Hamil- 
tonians (Sec. lIVF|) . We also use the example in Sec. lIV Al 
to discuss the convergence properties of the method 
(Sec. IIV Cj) . In addition to these one-dimensional exam- 
ples, we discuss the application of the algorithm to the 
two-dimensional Henon-Heiles system (Sec. IIV Gl) . 



II. POSITION-DEPENDENT MASS AND 
PT-SYMMETRIC SCHRODINGER EQUATIONS 

We want to replace m — > m(x) in one-dimensional 
quantum systems and, therefore, investigate various ways 
to order the operators in the kinetic-energy terms. A 
whole class of Hermitian Hamiltonians with position- 



2 



dependent mass is given by the von Roos Haxniltoniaxis 

m 

H = i (ri^prii^pm 7 + ih 7 p ih^prii") + V(x) (1) 

with 



a + /J + 7 = -1 



(2) 



and m = m(x) with the usual canonical operators p and 
x. The choice a = 7 = 0, /3 = — 1 leads to 



H= iplp + y(i), 

2 m 

while with a — — 1, /3 = 7 = one finds 

H = ifip 2 + P 2 i)+nx). 

4 V m my 



(3) 



(4) 



While there is no a priory reason to favor any of the 
Hermitian Hamiltonians, for example a specific choice 
of a, /3, and 7, non-Hcrmitian Hamiltonians are often re- 
jected for general reasons, since they can possess complex 
eigenvalues. However, in standard quantum mechanics 
closed systems are described by Hermitian Hamiltonians 
and real eigenvalue spectra. (Note that non-Hermitian 
Hamiltonians occur, e. g., in the approximate description 
of open systems with optical potentials.) Regarding the 
reality of the spectrum, non-Hermitian PT-symmetric 
Hamiltonians form a special class as discussed in [l2|, [l3| . 
Therefore, they are of special interest for possible exten- 
sions of standard quantum mechanics. Note, however, 
that these extensions often also involve a complex po- 
tential instead of a position-dependent mass. For the 
NH3 molecule the position-dependent reduced mass n{x) 
in [7| was derived within classical mechanics and then 
translated into quantum mechanics in the form 



H 



1 

2m" 



P 



V(x) 



(5) 



which is non-Hermitian, but convenient for computation. 
The advantage of Eq. (|5|) is that no first derivative of 
the wavefunction occurs. This simplifies the numerical 
solution with standard approaches. For example, the 
Numerov-Cooley method |T3, EH can be adopted as was 
done in Another possible non-Hermitian choice is 



H 



P 2 



1 
2m" 



(6) 



(Note if both the potential V and the mass m possess 
inversion symmetry, V(x) = V{~ x) and m(x) = m(—x), 
as is the case for ammonia, then the Hamiltonians in 
Eqs. (|5|) and j6]) arc P and T symmetric.) 

It is important to not only have an efficient and re- 
liable solver for one version of the Hamiltonian, but to 
be able to compare the results of different versions. Ide- 
ally, the results agree sufficiently well and thus the ques- 
tion of the proper form becomes practically inessential 



for the attempted approximation level. A strong varia- 
tion of the results with the chosen form of the Hamilto- 
nian, on the other hand, is a clear warning signal. Since 
non-Hcrmitian and (for example for V(x) 7^ V(—x)) in 
principle even PT non-symmetric Hamiltonians can be 
obtained for some versions of the Hamiltonians with a 
position-dependent mass, it is also of interest to have a 
solver that can handle these cases and that detects pos- 
sible non-vanishing imaginary components of the eigen- 
values. This is, in fact, also an important issue in the 
general context of PT-symmetric Hamiltonians (or other 
proposals for extensions of standard quantum mechanics) 
without position-dependent masses. Since in the major- 
ity of cases physically relevant Hamiltonians do not pos- 
sess analytical solutions, it is important to numerically 
check whether the eigenvalue spectrum is purely real or 
not. The algorithm presented below provides a promis- 
ing solution to this type of problems. Position-dependent 
masses are easily handled, even for different formulations 
of the kinetic-energy operator. Furthermore, a number of 
bound states are obtained simultaneously within a single 
calculation. Finally, the algorithm can equally well be 
applied to complex Hamiltonians and thus also to non- 
Hcrmitian Hamiltonians with either purely real or partly 
complex eigenvalues. 



III. THE MATRIX ALGORITHM 

The technique to investigate simple quantum systems 
reviewed below has been taught by one of us (U. W.) 
in the computational physics courses at Humboldt uni- 
versity since 2006. It has also recently been used for 
supersymmetric quantum mechanics in [f9j . The basis 
is the so-called SLAC derivative that was proposed for 
lattice fermions in (20|. While it had to be discarded 
for four dimensional quantum field theory due to incom- 
patibilities with ultraviolet renormalization plj no such 
problems seem to hamper the present quantum mechan- 
ical applications. 

We define a one-dimensional position space lattice con- 
sisting of an odd number of N = 2M + 1 points 



{-Ma,-(M- l)o,. 



{xi,x 2 , ■ ■ -,x N } = 
. , —a, 0, a, ... , Ma}. 



(7) 



that are equidistantly spaced with the lattice spacing a 
and lie symmetrically around the origin. Wave 'functions' 
in the Schrodinger representation of quantum mechanics 
are restricted to this space 



/ *(a:i) \ 
*(a*) 



(8) 



and thus the Hilbcrt space is approximated by C N . Intu- 
itively, a bound state that is centered around the origin, 
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can be represented well in this framework, if its size is 
much smaller than L = TV a such that the sites ±Ma are 
deeply in the classically forbidden region. Moreover a 
must be small enough to allow for a good resolution of 
structures in the wave function. For example, for a sim- 
ple harmonic oscillator of ma ss m and frequency u> these 
conditions amount to a <C yjfi/imuj) <C L. 

Linear operators must become finite matrices now. 
Functions of the position operator like the potential V(x) 
trivially translate into diagonal matrices 



V(x 2 )*(x 2 ) 
\ V(x N )V(x N ) J 



(9) 



To also implement a canonically conjugate momentum 
operator p that has to mimic the derivative, we need to 
impose boundary conditions which we take periodic with 
period L. Note that this also implies a periodic exten- 
sion of the diagonal elements of the position operator. 
This means that odd powers of x create jumps at odd- 
integer multiples of L/2, in particular at ±Ma, which 
will, however, be seen to cause no problems in our ap- 
plications here. The first idea that comes to mind now 
is to represent p 2 by a difference operator over three (or 
more) sites. This would clearly inflict leading discretiza- 
tion errors that are powers of a 2 . We found that a better 
precision is obtained by Fourier transforming, then using 
a diagonal operator for f> 2 analogous to the one in dis- 
crete position space and then transforming back. This 
results in a nonlocal matrix representation in position 
space that couples all sites given by 



(/(p))ifc = j; X] /(P) ex PM^ - x k)\- 



where the p-sum runs over the values 



(10) 



-M- 



,2-k 



{Pl,P2, 

x 2tt 

-{M-1) T , 



■ ,Pn} = 
,M— 



(11) 



and periodicity under shifts of p by ±27r/a holds. 

A multiplication of a wave function with the above ma- 
trix may be decomposed into two steps. First we compute 



JV 



*(p) = a exp[-ipxk]^(x k ) 



k=l 



L/2 



L/2 



dyexp[-ipy]^(y) 



(12) 



where the continuum limit a — > at fixed L is indicated 
(for a continuously defined In the second step we 
form 



(13) 



Here x can in any case assume continuous values and 
inside / obviously p is equivalent to —id/dx. The p-sum 
becomes infinite in the continuum limit. 

Returning to finite a and L the formula (setting j = 
i-k) 



(exp[iap]). tfc = -^^cxp[ip(a 



TV 
p 

(— iy sin(7ra/a) 



x k )} 



N sin(7r(a + ja)/L) 



(14) 



which holds for general a, allows for the explicit com- 
putation of matrix elements of powers of p. Comparing 
powers in a we find for example 







(P), 



for j = mod N 



ik 



(p 2 )* 



e iz^y pi SP 

f^(l - a 2 /L 2 ) for j = mod N 

2tt 2 (~1) J cos(7rj/JV) , 

TP- sin a (7rj/iV) elbt! 



(15) 



(16) 



Equivalently, we may of course also take matrix powers 

Of (p)jfc. 

To find the eigenvalues of the Schrodinger equation we 
now construct the matrix representation of H and diag- 
onalize this matrix to find eigenenergies E„ and eigen- 
functions \l/ n of the system. Clearly, the width L and the 
number of points TV should be sufficiently large to obtain 
converged results and only a certain number of low-lying 
energies E n can be expected to be reliable for any fixed 
TV and L. 

To diagonalize the matrix, we implement the matrix in 
MATLAB [22( and employ the routine eig which deter- 
mines eigenvalues and eigenvectors of a matrix (also for 
non-Hermitian matrices). A simple MATLAB code to 
construct the matrix representation of the Hamiltonian 
for the first application discussed in Sec. IIV Al (ammonia 
inversion) is explicitly given in the appendix. 

As will be discussed in Sec. IIV CI the matrix algo- 
rithm has a very fast convergence behavior such that 
around 100 points are usually sufficient. Therefore, also 
higher dimensional problems with for example about 
N x x TVj, = O(10 4 ) points in two dimensions are within 
reach on present-day computers. In that case we dis- 
cretize each direction as before and obtain a rectangle 
filled with sites 

fr = {xinVi*), ii = l,-.-,N x , i 2 = 1, . . . , N v (17) 

with the unique compound index 

I = h + (i 2 -l)N x = l,...,N x N y . (18) 

Operators are embedded in the tensor product state 
space. The position operators become for example 



(x) 



IK 



(y)iK = yi 2 $i 1 k 1 Si 2 k 2 - (19) 
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Similarly pj leads to a matrix 

(pDik = (p 2 )t 1 k 1 S, 2 k 2 , {vI)ik = Siiky (P 2 kfc 2 ( 20 ) 



with the matrices (p )ih taken from (|16p with the obvious 
substitution of a,L,N by the corresponding quantities 
referring to the respective direction. Note that in more 
than one dimension many vanishing matrix elements ap- 
pear in a typical Hamiltonian. One may thus consider to 
store it in a sparse matrix mode as a list of the nonvan- 
ishing elements rather than a full matrix. 

The two-dimensional Hcnon-Heiles system is treated 
as an example in Sec. HVGl and the MATLAB [H| code 
which generates the matrix representation of this two- 
dimensional problem is explicitly given in the appendix. 

The algorithm is not time critical for the one- 
dimensional problems discussed in this paper, which 
means we find well converged results for these systems 
in less than a second on a modern standard computer. 
For the two-dimensional Henon-Heiles system, around 5 
minutes are needed to find the eigenenergies and around 
60 minutes are required to find the eigenfunctions in ad- 
dition when using 101x101 points. Therefore, an im- 
plementation in programming languages like FORTRAN 
or C which would speed up the algorithm is a possi- 
ble future option, but was not needed for the exam- 
ples considered in this work. Furthermore, with only 
a single matrix diagonalization many bound states with 
eigenener gies E n are found simultaneously. Other al- 
gorithms jiol l23j which were previously used to solve 
the one-dimensional Schrodingcr equation with position- 
dependent masses consist of "guessing" an eigenenergy E 
and a subsequent test whether the resulting wave func- 
tion has the correct behavior to fulfill the Schrodingcr 
equation or not. 



IV. EXAMPLE APPLICATIONS OF THE 
MATRIX ALGORITHM 

A. The inversion motion of NH3 

To calculate the energy levels describing the inversion 
motion in ammonia with our algorithm, we exactly follow 
the procedure in m and first represent the data points 
(from table 1 in Q) for the double-well potential V(x) 
describing the inversion mode by an even polynomial of 
degree 20. Our coefficients are listed in the appendix. It 
might be noted that the potential obtained with DFT is 
in good agreement with the experimental value for the 
equilibrium geometry in p4| and the obtained barrier 
height of 2013.5 cm -1 agrees well with the one from em- 
pirical procedures (2018 ±10 cm -1 Q). In addition, the 
theoretical study in @ yields an effective barrier height 
of 2021 ± 20 cm^ 1 . The potential is therefore very useful 
for a one-dimensional study, though the method is nev- 
ertheless a little bit questionable, because the theoretical 
effective barrier in [8j is only obtained, if one takes into 
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FIG. 1: (Color online) Potential curve (solid red line, data 
from 0|) and wave functions (symbols, see legend in the 
graph) of the 6 energetically lowest states 0s to 2a. The wave 
functions have been shifted by constants such that their val- 
ues at \x\ — > 00 displays the eigenenergies. In the calculation, 
a width of L = 4 a.u. and N = 111 points were used. 



account the zero-point vibrational energies of the other 
vibrational degrees of freedom, which adds a value of 244 
± 14 cm -1 to the barrier height according to Q- 
The goal is now to solve the Schrodinger equation [37| , 



= 2fi(x) [V(x) - E] *(x), 

or equivalently, 



1 



p 2 + V(x) ) *(x) = EV(x) 



2fi(x) 

with the position-dependent reduced mass 



H(x) 



3mM 



3mx 



3m + M ri 



(21) 



(22) 



(23) 



where m is the mass of a hydrogen atom, M the one of 
the nitrogen atom, and r$ = 1.00410198 a.u. is given 
in Q as the N-H distance which minimizes the energy of 
the molecule in planar geometry. The result of the matrix 
diagonalization is illustrated in Fig. Q] and the obtained 
energies are listed in Table HI The states are named us- 
ing the symmetry label symmetric (s) or antisymmetric 
(a) together with an index n where n = stands for the 
energetically lowest state with symmetric or antisymmet- 
ric character. Thus 0s designates the symmetric ground 
state and we give the other energies relative to this level. 
A comparison of the energy levels found with the ma- 
trix algorithm to the energy levels in table U of shows 
perfect agreement. 

As has been mentioned in Sec.[TTl the obtained eigenen- 
ergies should ideally not change drastically, if one puts 
the position-dependent mass into the Schrodinger equa- 
tion in another form. Table [TT] shows the eigenenergies 
obtained with Eqs. © to ©. We find a maximal relative 
deviation of the eigenenergies listed in Table UJ of about 
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State 


/x=const [7] 


Ref. [7] 


This work 


Experiment 


Os 


0.00 


0.000 


0.000 


0.000 


Oa 


1.05 


0.837 


0.837 


0.793 [25] 


Is 


977.23 


931.72 


931.72 


932.43 [26] 


la 


1030.12 


968.67 


968.67 


968.12 [26] 


2s 


1651.69 


1596.76 


1596.76 


1598.47 [2JJ 


2a 


2011.44 


1885.33 


1885.33 


1882.18 [27] 


3s 


2558.75 


2389.14 


2389.15 


2384.17 [2JJ 


3a 


3142.66 


2902.99 


2902.99 


2895.61 [2JJ 



TABLE I: Inversion energy levels of NH3 (shifted by Eo B ) in 
cm -1 . The energies of our work were found with L = 4.0 a. u., 
N = 111 points and adopting the same ID potential as in Q. 



State 


Eq. © 


Eq. JH 


Eq. © 


Eq. © 


0s 


0.000 


0.000 


0.000 


0.000 


0a 


0.837 


0.833 


0.837 


0.837 


Is 


931.71 


932.01 


931.72 


931.72 


la 


968.64 


968.81 


968.67 


968.67 


2s 


1596.77 


1597.36 


1596.76 


1596.76 


2a 


1885.25 


1885.45 


1885.33 


1885.33 


3s 


2389.03 


2389.21 


2389.15 


2389.15 


3a 


2902.82 


2902.84 


2902.99 


2902.99 



TABLE II: Inversion energy levels of NH3 (shifted by Eo B ) in 
cm -1 for different implementations of the position-dependent 
reduced mass. In the calculations L = 4.0 a. u. and N — 111 
points were used. 



0.5% (0.005 cm -1 ) for the 0a state, indicating that the 
choice where to put the reduced mass has only a very mi- 
nor influence on the resulting energies within the given 
accuracy. 



State 


This work 


This work 


Experiment [25] 




/i=const Eq. (1241) 


fi(x) Eq. (HU) 




0s 











0a 


0.05 


0.05 


0.05 


Is 


793.8 


746.2 


745.6 


la 


798.3 


749.3 


749.15 


2s 


1419.7 


1368.4 


1359.0 


2a 


1513.7 


1432.0 


1429.0 


3s 


1912.6 


1836.4 


1830.0 


3a 


2238.4 


2106.4 


2106.6 



TABLE III: Inversion energy levels of ND3 (shifted by Eo a ) 
in cm -1 . The calculations were done with L = 4.0 a.u. and 
N = 111 points. 
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FIG. 2: (Color online) Convergence behavior with respect to 
a variation of the lattice spacing that is proportional to 1/N 
as we hold fixed the width L = 4.5 a. u. We show the ground 
state and the 7th as well as the 19th excited state. Relative 
errors of the energies l- E "^- Econv are plotted, where E conv is 
the 'converged' energy (average over the 10 largest N- values). 



B. Inversion energy levels of ND3 

In addition to NH3, we test our algorithm by finding 
the energy levels of ND3 which were not calculated in Q . 
This is also of physical interest, since it is a further check 
of the adopted position-dependent mass model. There- 
fore, we simply replace the mass of the hydrogen atom 
by the mass of deuterium m.D=2. 013553212712 amu [28[ 
and solve the Schrodinger equation ([5]). The obtained 
energies are listed in Table HTT1 It turns out that the in- 
version splitting can be reproduced and the energy levels 
for higher excited vibrational states have a relative er- 
ror ^~ Ecal °~' Ecxp < 0.7% in comparison to the experiment 
whereas the use of the constant reduced mass 

3mM ( 3m sin 2 B„\ , 
^=3^Tm( 1 + T^J (24) 

from HI with (3 e = 22° 13' (from Q) leads to higher 
energy values and, therefore, larger relative errors (up to 



6.6%). Even though the agreement with the experiment 
is very good, it still has to be noted that according to Q 
the barrier height depends on the zero-point vibrational 
energy of the other vibrational degrees of freedom, and 
this zero-point vibrational energy depends on the isotope. 



C. Convergence properties 

We discuss the convergence behavior of the matrix al- 
gorithm on the example of NH3 (Sec. IIV A[) . We use Eq. 
(|4]) as Hamiltonian and the number of points N is var- 
ied first for fixed L = 4.5 a. u. with the result shown 
in Fig. [2J In this way the dependence on the resolution 
given by the lattice spacing was tested. In Fig. [3] we 
have frozen the spacing to a = 4.5/151 a. u. and thus 
test stability when changing L. The nearly straight lines 
in the semilogarithmic plot indicate exponential conver- 
gence behavior with respect to both regulators. 

The convergence that we have demonstrated in our 



6 



LU 



LU 

I 



10" 



10" 



LU 10 



o 

LU 



10 



-10 



-15 



x\ 

• # x + 


• ground state 
x 7th excited state 
+ 1 9th excited state 


• ** T 

* X 4. 
• X + 

• X 

• + 

• X xkiMtau 





100 120 140 
Number of points N 



160 



FIG. 3: (Color online) as Fig. [2] but now a = 4.5/151 a. u. is 
fixed and thus L varies proportional to the number of points 
N. 



first test is actually physically plausible. As indicated be- 
fore the only exponentially small sensitivity with respect 
to (large) L is due to the smallness of the wave functions 
for bound states deep in the classically forbidden region. 
The same wave functions in momentum space will also 
have an exponential fall-off at large momenta so that a 
dual argument holds in momentum space. We may argue 
here by analogy with the sampling theorem [3uT ]. A con- 
tinuous time signal that contains no frequencies beyond 
the Nyquist frequency u = tt/t can be exactly recon- 
structed from sampling it at discrete times separated by r 
(CD player). In analogy, if the support of our bound state 
wave functions would be exactly contained in the interval 
[— 7r/a, 7r/a] then the exact wave function for continuous 
x could be reconstructed form the Fourier components 
(discrete and finite in number for finite L). Then there 
obviously is an exact correspondence between derivatives 
d/dx and factors ip. In reality the boundstates do not 
have compact momentum support, but the deviation is 
only caused by the exponentially small talcs for small 
enough a. 



D. Morse potential 

The discrete position (momentum) space introduced in 
the matrix algorithm is symmetric with respect to x = 
(p = 0) as is also the double- well problem discussed in 
the previous section. Nevertheless, general Hamiltonians 
can be non-symmetric and may also have a partially con- 
tinuous spectrum (in the infinite volume). Therefore we 
want to investigate, if bound states can still be found ac- 
curately and efficiently using the matrix algorithm also 
in this more general case. A popular example for a non- 
symmetric potential with bound and continuum states is 
the Morse potential 



V(R) = D e (l 



(25) 



Often the vibration of a diatomic molecule can be well- 
described by this potential. The energy spectrum for the 
bound states of this potential is known analytically HJ , 



D e 



(26) 



with 



where fi is (again) the reduced 



_a_ / 2D B 
2tt y n 

mass. We investigate a Morse potential with parameters 
D e = fx = 1 and a = 0.24, R e = —35 where 6 vibrational 
bound states exist. Table LTVl shows the energies obtained 
by the matrix algorithm with N = 111 points and L = 90 
together with the exact results from Eq. (|26|) . We still 
find excellent agreement except for some visible discrep- 
ancy for the highest lying 6th vibrational state. This can 
be understood, since states close to the continuum are 
more extended in space. UN — 201 points, L = 140, 
and R e = —60 are used, all eigenenergies obtained by 
the matrix algorithm are identical to the analytical re- 
sults within the accuracy given in Table IIVI 

In addition to the bound states, we obtain discretized 
continuum states on which the periodic boundary con- 
ditions of the method arc imprinted. To check whether 
this discretized continuum can approximately represent 
the true continuum, we examined as an example the re- 
lation 



(0|x 2 |0) « <0|x|n> (n|x|0) 



(27) 



n=0 



using the discretized states \n). Using N = 301 points, 
L = 140, and R e = —60 the relative error 



^max) 



<0|x 2 |0>- £ (0|x|n) (n|x|0) 

71 = 



(0|x 2 |0) 



(28) 



was calculated. With a restriction of the summation to 
bound states only we find e(n max = 5) w 10~ 6 . When 
we further increase n max the error decreases monotoni- 
cally until machine precision is reached (e(n max ) < 10~ 14 
for rt max > 116). Thus the completeness relation of the 
eigenstates is numerically fulfilled, if also the discretized 
continuum states are considered. We conclude that any 
upcoming non-symmetric problem that contains a con- 
tinuum spectrum should also be well treatable by the 
matrix algorithm, if discretized continuum states (with 
periodic boundary conditions) are sufficient as is the case 
in the here considered example of bound to continuum 
transitions. 



E. Harmonic oscillator with a position-dependent 
mass 

To compare the present matrix algorithm to the re- 
cently published traditional second-order Hartrcc shoot- 
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n 


Matrix algorithm 


Exact 





0.1625056275 


0.1625056275 


1 


0.4443168825 


0.4443168825 


2 


0.6685281374 


0.6685281374 


3 


0.8351393923 


0.8351393924 


4 


0.9441506473 


0.9441506474 


5 


0.9955620565 


0.9955619023 



TABLE IV: Eigenenergies of the Morse potential ((25j ob- 
tained with TV = 111 points and L = 90 for the model param- 
eters D e = h = fi = 1 and a = 0.24, R e — —35. The exact 
values were obtained using Eq. (|26|) . 



ing method [T^, we implemented the two model Hamil 
tomans 

1 . 1 . 1 .o 



Hi = ^ P 



2- 1+^ P+ 2 X 



and 



w 1 - f 1+i 



x 2 . 



(29) 



(30) 



When we compare our eigenenergies obtained by the ma- 
trix algorithm with TV = 201 points and L = 20 (which 
is well within the converged regime), we reproduce ta- 
ble 1 (for Hi) and table 3 (for H 2 ) in [l^ to all digits. 
Considering the convergence properties, qualitatively the 
same exponential behavior is found as shown in Figs. [5] 
and One should note that the second-order shoot- 
ing method in 

[3 only converges quadratically and, 
therefore, much larger numbers of evenly spaced points 
(TV = 640, 960, 1440 and 2160) were needed to obtain the 
results. 




FIG. 4: (Color online) Top: Relative Error of the real part of 
the eigenvalues of Eqs. (|31|) [x] and (|33|) [+]. The values were 
obtained with the matrix algorithm using TV = 101 points and 
L = 25. The exact values E n ex are given by Eqs. (|32[) and 

HMD- 

Bottom: Absolute Error of the imaginary part of the eigen- 
values (same L and TV). 



has according to [12| the complex eigenvalues 



F. Non-Hermitian PT-symmetric and 
non-symmetric cases 

As mentioned in Sec. UH the PT-symmetric Hamilto- 
nians examined in Sec. [XT] and IIV Al are rather special. 
Therefore, we will briefly discuss the more usual form 
of a PT-symmetric Hamiltonian with a constant mass 
but a complex potential. We want to demonstrate that 
the matrix algorithm is also applicable for such types of 
problems, since complex matrices can also be treated. 
Therefore, we implemented the model Hamiltonian 

H = p 2 + x 2 + i x. (31) 

According to [l2[, this PT-symmetric Hamiltonian has 
the real eigenenergies 

E n =2n+^. (32) 

In contrast, the non-Hcrmitian and not PT-symmetric 
Hamiltonian 

H = p 2 + x 2 + i x - x (33) 



E n =2n + l + -i (34) 

which can no longer be interpreted as eigenenergies. The 
relative error of the real part as well as the absolute er- 
ror of the imaginary part of the eigenvalues obtained with 
the matrix algorithm with TV = 101 points and L = 25 
is shown in Fig. |4l The eigenvalues of the Hamiltonian 
in Eq. pip now contain a small imaginary part (also 
see Fig. 2]) due to roundoff errors which was not the 
case for the PT-symmetric Hamiltonians in Eqs. (0 and 
©. The reason is that the matrix representations of 
Eqs. ([5]) and (|6]) are real, while the matrix representa- 
tion of Eq. (|3"Tj) is complex. Nevertheless, the very good 
agreement of our results with Eqs. (|32|) and (|34|) (rela- 
tive error below 10~ 12 for n < 45) shows that our matrix 
algorithm can also be very helpful to find eigenenergies of 
PT-symmetric Hamiltonians. In critical cases one could 
in principle consider increasing the numerical precision 
(cither in steps from single to quadruple precision or, 
even better, digit-by-digit) to probe further whether the 
eigenvalue spectrum is purely real or not. 
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-10 -10 



FIG. 5: (Color online) Potential energy surface (filled white) 
and obtained (discrete) ground state wavefunction (filled blue, 
scaled arbitrary) of the Henon-Heiles system in Eq. (|35[) with 
A = 1/V80. A L x x L y = 20 x 20 grid with N x x N y = 61 X 61 
points was used. 



G. The two-dimensional Henon-Heiles system 

To demonstrate the performance of the matrix algo- 
rithm in two dimensions we investigate the Henon-Heiles 
system [32[ 



This model potential is frequently used as a benchmark 
for numerical methods [33l - l35j , although the potential in 
Eq. (f3"5|) is not bounded from below and thus it does 
not support true bound states but only metastable ones 
that decay by tunneling through the barriers as is dis- 
cussed in, e.g., Refs. |35|, |36[. We use A = l/y/80 to 
compare our lowest 36 eigenenergies with the ones ob- 
tained in [H]. The MATLAB [H program which gener- 
ates the matrix representation of Eq. Q35p is given in the 
appendix. When using a L x x L y = 20 x 20 grid with 
N x x N y = 61 x 61 points, we reproduce all 36 eigenen- 
ergies given in (35j within the given accuracy. The ob- 
tained ground state wavefunction is shown together with 
the potential in Fig. [5] When increasing the number of 
grid points to N x x N y = 81 x 81 and N x xN y = 101 x 101 
(L x x L y = 20 x 20 held constant) or changing the grid 



size to L, 



Ly = 



18 x 18 and L x x L y 



22 x 22 



(a x = a y = 20/61 approximately held constant), the 
obtained energies are stable within an accuracy of at 
least 12 significant digits (38|. Thus the matrix algo- 
rithm is also very suitable for finding bound states of 
two-dimensional Hamiltonians. 



V. SUMMARY 

The matrix algorithm presented in this paper is easy 
to implement, flexible, and shows exponential conver- 
gence (with respect to the number of g rid points N and 
width L). While other algorithms [Io[, HH usually con- 
sist of guessing an initial energy, this algorithm repre- 
sents a much more direct method, since one finds a set 
of eigenenergies by just a single matrix diagonalization. 
To demonstrate the performance of the algorithm it was 
applied to a one-dimensional model describing the in- 
version motion of NH3 and ND3. Furthermore, differ- 
ent forms for the position-dependent mass Hamiltonian 
were discussed. A great advantage is the flexibility of 
the present algorithm to handle their different possible 
forms. Especially, it is possible to avoid non-Hermitian 
model Hamiltonians (as used for, e.g., describing inver- 
sion of NH 3 ) that are often only adopted, because the 
symmetrized Hcrmitian form of the Hamiltonian leads to 
a Schrodinger equation that cannot be solved with many 
of the standard numerical algorithms. 

In the case of ammonia it turns out that the eigenen- 
ergies do not strongly depend on the choice where one 
puts the mass into the Hamiltonian. This justifies the 
underlying classical derivation of the model Hamiltonian 
with a position-dependent mass. Also the eigenenergies 
of a Morse potential were determined. The results indi- 
cate that the algorithm is applicable to find bound states 
of general Hamiltonians that may contain asymmetric 
potentials and include a continuous spectrum. The al- 
gorithm was also successfully applied to non-Hermitian 
Hamiltonians both with or without PT symmetry. This 
allows to numerically check whether the eigenvalue spec- 
trum of such an Hamiltonian is purely real or not and 
demonstrates that the algorithm can handle Schrodinger 
equations with complex eigenvalue spectra. Finally, it 
was demonstrated with the aid of a two-dimensional 
problem that the algorithm is straightforwardly applied 
to higher-dimensional problems, only limited by the in- 
creasing numerical efforts due to the exploding number 
of grid points needed. 

For very large problems or to go beyond two dimen- 
sions it may be useful to write a routine that applies H 
to a vector 'F and renounce at storing any matrix at all. 
Such a routine would apply FFT twice for each direc- 
tion in every call with the cost scaling like iV^lnA^. If 
only some low- lying energies are requested, the use of a 
Lanczos- or Arnoldi-type algorithm would be a natural 
choice for enhancing the efficiency Within Matlab the 
eigenvalue finder eigs would allow for such an approach 
and f f t is available as well. 



Acknowledgments 



Wc thank Dr. Oliver Bar for helpful discussions. 



9 



Appendix A: Matlab code fragments 

The following short MATLAB [22| function generates 
the lattice and the real antisymmetric matrix if> for given 
integer M and width L that will enter into the construc- 
tion of Hamiltonians: 

function [x,ip_op] = gen_mommatrix(L,M) 

'/ generate the lattice x (vector) 

'/ and i times p-operator (real, antisymmetric) 

N=2*M+1; 

a=L/N; '/ spacing 

x=a*(-M:M); '/, x-values 

ip_op=zeros(N,N) ; '/ matrix p 
cl=pi/L; 
c2=pi*(N+l)/N; 
for i=l:N-l 
for k=i+l:N 

ip_op(i,k)=cl/sin(c2*(i-k)) ; 

ip_op(k, i)=-ip_op(i ,k) ; 

end 
end 

For the ammonia inversion problem the following con- 
stants arc required: 



1.1565093681631e+02 
-3.2323821164423e+02 

5.4331165379878e+02 
-5.0630533518111e+02 

2.0128292638493e+02 ]; 



The Hamiltonian is constructed by the following sequence 
of steps (for L = 4, N = 111): 

[x,ip_op] = gen_mommatrix(4,55) ; 
'/, x-dependent mass 

mux=3*m*M/ (3*m+M) +3*m*x . "2 . / (r0~2-x . ~2) ; 
mux=mux*amu; °/ conversion amu -> a.u. 
'/, kinetic part 
psq=-ip_op~2 ; 

H=-0 . 5*ip_op*diag(l . /mux) *ip_op; °/Eq. (3) 

y„H=0.25*(diag(l./mux)*psq+psq*diag(l./mux)) ;'/Eq. (4) 
y„H=0.5*diag(l./mux)*psq; '/Eq. (5) 

'/„H=0 . 5*psq*diag(l . /mux) ; °/Eq. (6) 

y, add potential part 

Vpot=polyval (f liplr (K) , (x*aBohr) . ~2) ; '/ evaluate 

% polynomial 

H=H+diag(Vpot) ; 

Eigcnencrgies and cigenfunctions can now be found using: 



'/.constants from \protect\vrule widthOpt\protect\href -[htt^W?t/p£uBi8sElllg€if 6H/cuu/}-[http : //physics .nist . gov/ cuu/} 
Hartreecm=219474. 63137;'/! a.u. (Hartree energy) 



aBohr=0. 52917721092 



"/constants from ref 
m=l. 007825035; 
ym=2. 013553212712; 
M=14. 003074; 
amu=1822.888; 
r0=l . 00410198/aBohr ; '/N-H-distance in 
y, Potential fit parameters: 
y„V=sum_{i=0}~{i=10} K(i+1) x_~{2i} 
'/, [x in Angstrom, V in a.u.] 
K=[ 

-1 . 2760373471398e-01 

4.7973549262032e-01 
-4 . 4967805753691e-01 

3.4048981035460e+00 
-2 . 5268066877745e+01 



"/.in 1/cm 

Vol a.u. (Bohr radius) 
°/in Angstrom 
Aquino et al . 
'/hydrogen mass in amu 
'/deuterium mass in amu (from nist) 
'/nitrogen mass in amu 
°/l amu in a.u. 



The following MATLAB [22| program generates the ma- 
trix representation for the two-dimensional Hcnon-Hcilcs 
system fl35]) using N x = N y = N = 2AI + 1 points with 



L, 



L: 



[x,ip_op] = gen_mommatrix(L,M) ; 
psq=-ip_op~2 ; 

del=eye(2*M+l) ; '/ unit matrix 
H=0.5*(kron(psq,del)+kron(del,psq)+. . . 

kron(diag(x. ~2) ,del)+kron(del,diag(x . ~2) ) )+ . . . 

(l/sqrt(80))*(kron(diag(x. ~2) ,diag(x))-. . . 

kron(del,diag(x. ~3))/3) ; 

We have used the MATLAB function kron that builds 
the tensor (Kronecker) product of two matrices. The 
generalization to discretizations that are anisotropic in x 
and y and also to more than two dimensions is obvious, 
and eigenvalues of H are found as before. 
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